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obtained. We compare the mean-field approximation and the continuum approximation 
to the exact solutions and we discuss their regime of validity. 
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I. INTRODUCTION 

A standard way of studying the kinetics of reaction-diffusion systems involves deriva- 
tion of an infinite set of moment equations for the state probability. The solution of such 
equations usually poses great mathematical difficulties; however, the simple topology of 
the one-dimensional lattice frequently allows derivation of exact solutions These 
results may be used to test the validity of various approximations. Neglecting spatial 
fluctuations in concentration and occupation number space leads to classical macroscopic 
rate equations. This type of approximation, which implicitly assumes that each particle 
interacts with the system as a whole through a mean field (MF) , should improve as the 
number of interacting neighbors grows, i.e., with increasing dimensionality of the lattice. 
Although this approach may seem too rough for one-dimensional systems, it works well 
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for early time evolution and in other situations where the distribution of particles is 
nearly free of correlations. 

Exactly solvable reaction-diffusion models consist largely of single species reactions 
in one dimension, e.g., variations of the coalescence process, A + A —>■ A + S [f7|-[T2| 



and the annihilation process A + A — > S + S |T0,O,O-21], where A and S denote 



empty and occupied sites, respectively. These simple reactions display a wide range 
of behavior characteristic of non-equilibrium kinetics, such as self organization, pattern 
formation, and kinetic phase transitions. Interval methods have provided many exact 
solutions for one- dimensional coalescence and annihilation models. The method of empty 
intervals, applicable to coalescence models, requires solution of an infinite hierarchy of 
differential difference equations for the probabilities En of finding n consecutive lattice 
sites simultaneously empty |7|-p|, p!2| , p2| -p^ . For annihilation models, the method of parity 
intervals similarly requires determination of the probability of n consecutive lattice 
sites containing an even number of particles p5|-P^. 

The infinite hierarchy of equations for the interval probabilities can be used as a 
starting point for the derivation of MF and continuum approximations. Most results ob- 
tained from interval methods rely on a continuum approximation of the spatially discrete 
equations; the resulting PDEs allow for straightforward solution. However, the pertinent 
set of differential difference equations for the interval probabilities admits analytical solu- 
tions ||2^, providing exact results which may differ from the continuum limit. In general, 
we expect that this discrepancy becomes important for high concentrations and short 
time scales (compared to the typical reaction time), since on-lattice local concentration 
perturbations travel at finite speed, in contrast to propagation in a spatial continuum. 

In this work, we use the interval methods described above to study the on-lattice 
kinetics of coagulation and annihilation reactions including particle creation steps. Ex- 
plicit expressions for the lattice coverage and the concentration are derived and compared 
to those yielded by the MF and the continuum approximation. The behavior of coales- 
cence and annihilation changes qualitatively if the reaction schemes incorporate particle 
sources in the form of input (external source) or back reactions (internal source). In the 
absence of sources, the concentration displays an anomalous t~^/^-decay to an empty. 



steady state ||T8|J19|] . Particle sources give rise to steady states with nonzero concentra- 
tion and change the transient behavior PJT5[]. In most cases, the concentration may be 
thought of as an order parameter in a "phase transition" between the empty and the 
active steady state, which is controlled by the effective rate of particle creation, h. In 
analogy to the theory of critical phenomena, one can regard h = as a transition point 
and characterize the steady state and the relaxation time of the system near the critical 
point through a set of static and dynamical exponents p9| . We shall see that our method 
allows the characterization of the system not only near the phase transition, but also far 
beyond criticality. 

This work is organized as follows. In section |I|, we define the general form of the 
coalescence model and solve for the cases in which the back reaction or a homogeneous 
particle input is present. The validity of MF-like rate equations and the continuum ap- 
proximation is analyzed vis a vis the exact solution. In section we do the same for the 
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annihilation model with different Icinds of particle input. Similarities with Glauber-type 
spin models and the relevance of these models to physical systems are briefly discussed. 



We also treat an annihilation model with symmetric birth |^ , both with immobile and 
diffusing reactants. We conclude in section with a summary of our results and possible 
extensions of this work. 



II. COALESCENCE MODELS 



All reactions in sections |l| and |T| are defined on an infinite ID lattice with spacing 



a. The basic coalescence model involves particles moving randomly and asynchronously 
to nearest neighbor sites with hopping rate 2D/a^\ in an extended model, particles give 
birth, i.e., a particle at a given site generates offspring at an empty adjacent site at rate 
v/a. One can also include a homogeneous source: particles are injected into the lattice 
at rate Ra. In all particle disappears whenever it lands on another. Let En{t) 

be the probability that a randomly chosen segment of n consecutive sites contains no 
particles. By noting the changes in En during a small time interval. At, and taking the 
continuous time limit At 0, we derive the Master equation 0: 

— TT ^ -^{En+i — '^En + En^i) + -{En+i — E^) — RnaE^. (1) 
at a 

The first term on the rhs represents the effect of the net particle flow into and out of 
an empty interval, whereas the second and third terms describe the effect of the back 
reaction (cooperative particle birth) and a homogeneous particle input, respectively. By 
comparing (|l]) for = 1 to the changes in Ei during a time interval At, one arrives at 
the boundary condition -Eo(^) = 1- Also, for non-empty lattices, Eoo{t) = 0. The case 
D > 0,v = R = has already been solved in a previous paper [^]. In that case, the 
reaction displays universal, anomalously slow t"^/^-decay to an empty state, as opposed 
to the MF t~^ asymptotic decay. This result holds both on-lattice and in the continuum 



limit p8[, and is reminiscent of persistent transient fluctuations induced by the dynamic 
self-ordering of the system 0]. 

In the sequel, we focus on the cases D > 0,v > 0, R = (reversible coalescence 
A + A ^ A + S, subsection |ll7\| ) and D > 0,v = 0, R > (coalescence with input. 



subsection [II B| ). We use the initial condition 

E„(0) = (l-po)^ (2) 

This corresponds to a random homogeneous distribution of particles characterized by a 
global coverage po and a concentration (number of particles per unit length) cq = po/a. 
Our primary goal is to compute the time evolution of the global coverage p(t) = 1 — Ei{t) 
and the associated concentration c(t) = p{t)/a. 
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A. Reversible Coalescence 



The existence of a kinetic phase transition for this system in the continuum hmit is 
well known Recently, Lin has shown that the phase transition is not an artifact of the 
continuum limit, since asymptotically this approximation only gives rise to quantitative 



corrections |]3T|. To do so, he solved the corresponding empty interval hierarchy by 
mapping it into a spin system with Glauber dynamics. However, it is possible to solve the 
hierarchy more straightforwardly by means of a simple Laplace transform (LT) method, 
as we show below. 

The boundary value problem (EVP) of interest is ||^ 

dEn 

—, — = En+l — 2En + En-1 + h {En+1 — En) , (3) 
CLT 

where r = 2Dt/a^ is a dimensionless time and h = va/2D is the relative feed rate; 
the boundary conditions read Eq{t) = 1 and E^{t) = 0, and the random homogeneous 
initial condition is that given by First, we solve for the steady state En,s- Setting 
dEn 

= on the Ihs of Eq. (0) , we get 

dr 

(1 + h)En+l,s - (2 + h)En,s + En-l,s = 0. (4) 

This difference equation has the solution 

En,s = (1 + h)-", (5) 
whence the steady state concentration is obtained: 

c. = ^ = l^ = — (6) 
a a 2D + va 

It has been shown that this is a true equilibrium steady state characterized by the 
presence of detailed balance and a maximum of the entropy 0]. The spatial coherence 
induced by the coalescence reaction is eventually destroyed by the back reaction. 

Eq. may be solved by standard techniques, for example, by taking the LT with 
respect to r, fitting a power ansatz to the resulting difference equation, and finally 
inverting the LT of the solution. With the boundary data -Eo(^) = 1 and Eoo{t) = 0, the 
LT of Ei(t) is given by 



where /3 = hpo — pl/{l — po)- Using the faltung theorem for the LT this can be 
inverted to yield 
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1 - Po - 



where /«(■) is the n-th order modified Bessel function. If the lattice is completely filled 
initially (po 1), the difference equations for En{s) are homogeneous and the second 
term on the rhs of (H) vanishes. 

Using the s oo expansion of Ei facilitates investigation of the short time behavior: 



EAs) 



1 - Po I3{1 - Po) 



+ 0{s' 



By virtue of a Tauberian theorem p3| , |3^ , we obtain 

p{r) = Po + \hpQ{l - Po) - Pol ^ + 0{t^) 



r ^0, 



(9) 



(10) 



for the short-time site occupancy. 

In the opposite limit, the expansion of Ei{s) about s = provides no additional 
information about the long time solution; only the steady state solution is recovered. 
Determining the long time asymptotics requires a bit more work. We write Ei{t) as 



2 r'^ I rlr' 
Ei(r) = l-p,-- / e--^/i(7r') — 
7 Jt t 

7 Jt T 



-/3t 



111 



with the notation 

a = 2 + h 7 = 2Vl + h ps 



l + h 7 



7 JO r' 



:i2) 



All integrals in Eiij) converge, since a — /3 > 7. Moreover, [/ is a known LT: 

2 



U 



r 



(13) 



Based on the value of the initial coverage, po, there are three cases for U: 



U 



1 - Po for Po > Pc, 
for Po = Pc, 
for Po < Pc, 



/h+l 



1 



(fe+l){l-po) 



where p^ = 1 — ^/T^^])^ is a critical initial coverage. One has pc < Ps, and pc 



(14) 



2Ps as 
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For long times, the asymptotic form of Ii{z) may be used to approximate the 
remaining integrals in the expression for -E'i(r). We then have the following asymptotic 
formula for the coverage: 

(Q-/3-7){a-7) tor po > Pc 

-^r-'/'e-^"^-^^- for po = Pc (15) 

- (i - Po - ih+m-po) ) e"^" Po < pc- 

Considering then the (dimensionless) relaxation time tr, defined as 

= - jH^^"^ln|p(r) -p,|, (16) 

we see that the system undergoes a second-order dynamical phase transition, from a 
slowly relaxing phase — with Tpt dependent on po — for po < Pc, to a phase with 
constant r^j, for po > Pc- It can be readily shown that the first derivative of tr at po = Pc 
is continuous, while the second derivative displays a finite jump determined by the value 
of h. 



MF and continuum approach 

We now examine the simplest type of MF approximation for reversible coalescence. 
In this approach, the kinetics of p may be obtained either by means of simple heuristic 
arguments or directly from the empty interval hierarchy. Let P„ be the probability of n 
consecutive sites being simultaneously occupied. Making use of the identities Ei = 1 — Pi 
and = 1 — 2 Pi + P2 in the evolution equation for Ei, we have 

^ = -il + h)P, + hP,. (17) 
dr 

Neglecting spatial correlations amounts to setting P2 = Pi in this equation. Taking into 
account that Pi is identical with the mean site occupation p, we get: 

^ = -^l + h)p'' + hp. (18) 
ar 

The physically acceptable steady state of this equation is given by 

^ P. ^ ^. (19) 

The time dependent solution matching this steady state is 

^ [{h/po)-l-h]e-'^- + l + h- ^ ^ 
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For short times, the Taylor expansion of the exponential function yields (]T0|); as expected, 
the MF approach is a good approximation at early times. In the opposite, long-time limit, 
an exponential relaxation towards the equilibrium steady state occurs: 



MF 

p - Ps (X Ps 



1-^ 



e"'^^, r ^ oo. (21) 



Po. 

Thus, in the MF approach, is given by 

TR = h-\ (22) 

This result is universal in the sense that the relaxation time does not "remember" the 
initial condition po, which is only contained in the prefactor 1 — ps I po- 
ll the initial coverage is low enough, the relaxation towards equilibrium is basically 
given by the time taken to fill large gaps through particle birth from the edge sites IQ]. 
Therefore, becomes initial-condition dependent in this regime, in contrast to the MF 
result (c.f. Fig. |l]). Moreover, tr becomes arbitrarily large in the limit po ~^ 0. On the 
other hand, if po > Pc, tr is MF-like only for large h. Note, however, that the relaxation 
is not purely exponential, since additional powers of r appear in the asymptotic form of 



the exact solution of Eq. (|T^). 

The exact results can also be compared with the continuum approximation, where 
one lets the lattice constant a shrink to zero To this end, we rewrite (|l^) in terms 
of concentration, using Cg = (1/a) h/{l + h) and 2D/v'^ as the units of concentration and 
time, respectively: 

c = co+ "°^\~"°^ t + 0(t^), (23) 
h 

where c and t are dimensionless now. On the other hand, rescaling c and t as above and 
taking the limit h ^ for fixed R and D in the long time form of the solution, Eq. (|15]), 
yields 

1 + 27r-V2[i _ (1 _ 2co)-2]t-3/2e-*/4 for cq > 1/2, 
1 - 7r-i/2ri/2e-</4 for Co = 1/2, (24) 

1 - (1 - 2co)e-^«(^-^o)* for < cq < 1/2. 

An alternative way to perform the continuum approximation consists of replacing the 
set of difference equations (^ by the PDE 

^^.2D^^.f, (25) 

ot ox^ ox 

for a spatially continuous function E{x, t), the boundary conditions now being E{0, t) = 1 
and E{oo,t) = 0. This equation can be solved for E{x,t) with a suitable exponential 
ansatz 0; the global concentration c{t) is then obtained by deriving with respect to 
space: 
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c{t) 



dE{x,t) 
dx 



(26) 



x=0 



Eventually, one can rescale c and t in dimensionless units as done above to obtain 



c{t) = 1 - -erfc i-Vi 



Co 



-Co(l-Co)t 



erfc 



- Co 



, (27) 



where 6(-) is the Heaviside step function. As expected, the long time development of 
this expression is identical to (p^. In contrast, for short times one has 



c(t) = Co + 



,(1/2 -co)^ 



.2A 



Vt + (1/2 - co)co(l - co)t + Oif'^). (28) 



Thus, for early times there is a qualitative discrepancy with respect to the on-lattice 
case, even as co — 0. The evolution is initially faster in a continuum, due to the infinite 
speed of propagation of local perturbations. 



B. Coalescence with input 



As a second example, we consider the effect of a homogeneous external input S A, 
at rate Ra, on the simple coalescence reaction. This model is described by the set of 
equations [0: 



dr 



Er,.i - {2 + hn)Er, + En+i, 



(29) 



with the boundary conditions Eq = 1 and E^o = 0, where we used the dimensionless 
time r = 2Dt/a? and relative feed rate h = Ra? /2D. 

We now compute the exact steady state solution by setting the Ihs of (p9[) equal to 
zero. The resulting set of difference equations can be compared to recursion relations 



satisfied by the Bessel functions [35 



JU-1{Z) - ^ Jy{z) + Ju+l{z) 



0. 



(30) 



The boundary condition at the origin determines the appropriate normalization factor 
for the steady state solution, En^s- The steady state coverage is then given by 



Ps 



2/1-1+1 



J2^-i (2/1-1) 

For /i <^ 1 we can use the Taylor expansion of Jx+i\x{x) around x to obtain 

1 d 



(31) 



Ps 



J2/i-i(2/i 1) dx 



(32) 



x=2/i-l 
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We compute the derivative of the Bessel function by using the asymptotic form p5| 

J,(z/ + zz/^/=^) = 2^/V-i/3Ai(-2i/=^2) + 0(z/-^), z/^oo, (33) 

(where Ai(-) denotes the Airy function) and obtain an exphcit expression for the small 
h limit of the stationary coverage: 

h-.^ (34) 

In the opposite limit, h ^ oo, we have the following result: 

Ps = 1 — h^^, h ^ oo. (35) 

Let us now examine the early time kinetics of our model. In general, it is difficult to 
obtain a closed solution of the hierarchy ( pQ] ) because of the nonconstant coefficient of 
En on the rhs. Under the simplifying assumption of an initially full lattice, we take the 
LT of (g|) and get 

En-i -{2 + S + hn)En + En+i = 0, (36) 

with Eq = 1/s and E^o = 0. Exploiting again the analogy of Eqs. (p6D with the recursion 
relations for the Bessel functions and using the boundary conditions leads to 

P _ 1 ^(2+s)fe-i+n(2/^~^) , - 

. J(2+.).-i(2/.-i) • ^'^^ 

For s 0, one recovers (^) by making use of the theorem lims^o sEn = -E„(oo). Since 
the conjugate time variable s appears in the index of the Bessel functions, the exact 
inversion of (P7| ) is quite involved. Instead, it is better to work with asymptotic formula. 
For n = 1, we can express Ei as a continued fraction |35|| : 



E,{s) = (-) ^ . (38) 

Vs/ s + 2 + h i J 

s + 2 + 2h 

s + 2 + 3h- ... 

We can now expand the right-hand side in powers of by using the Stieltjes method 
for J-fractions |^| and find 

^.w4-^.^.-(i). 

Inverting this expression term by term, we obtain the early time behavior of the empty- 
site coverage (c.f. Fig. H): 

p^r) = l-T+(l + ^] r2-i^±^!±ir3 + 0(r^), r ^ 0. (40) 
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The continued fraction representation of the LT is not suitable for obtaining the long 
time asymptotics. Instead, we make an exponential ansatz directly in Eqs. (^9t): 



= ^„,. + a„e-2^^ (41) 
The calculation of follows the guidelines laid by Racz for annihilation with input 



T5[| and paragraph |III A| ) . Substituting Eqs. pl|) in the set of differential- difference 



[see 

equations (|29|) , one finds that the coefficients a„ obey the same set of equations as En, 
except that 1 +nh/2 is now replaced by 1 — X + nh/2. Using the boundary condition at 
n = 0, we find the equation for the possible values of A: 

J(i-A)2h-i {2h~') = 0. (42) 

All the solutions Xi{h) are positive, since the zeros of the Bessel functions satisfy the 
inequality u < meaning that all homogeneous perturbations around the steady 

state decay exponentially with time. The small h limit is carried out using the asymptotic 
form (H). Setting z/ = 2(1 - X)h-^ and z = v'^l'^jiX - X) yields 

kx{-2}r''I^XI{\-Xfl^)=^. (43) 

From this, a cubic equation for A can be extracted, 

A?+A-P = 0; V=^^h\ (44) 

where are the zeros of the Airy function. Since the discriminant A = (p/3)'^ + (p/2)^ 
is positive, this equation has two unphysical complex solutions and one real solution, the 
latter of which is given by Cardan's formula: 

A,=M + w; u = (p/2 + VA)^/^ t; = - — . (45) 

u 

Expanding this expression for small /i, we find the smallest eigenvalue 

Ai = ^ h^l\ h 0, (46) 

from which we obtain 

= 2Xi = \ai\h^/^ + 0{h^/^). (47) 

The corresponding MP equation for p(r) may be obtained following the procedure of the 
preceding subsection. We find 

^ = V + Ml-P). (48) 
Integrating and using the initial condition, we obtain 
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where r] = y/W^^Ah, p± = {—h±ri)/2 and a = (po — P+)/(po — P-)- The predicted long 
time coverage is p'^^{oo) = pf^ = p+. For small h, the steady state coverage is 

pf^ = h'/' + 0{h), (50) 

whereas for large h we get 

pf^ = l -/,-! + 0(/^-2)_ (51) 

We now examine the transient behavior of the MF approximation. Taking po = 1, i.e., 
an initially full lattice, the early time kinetics is 

(52) 

while the long-time behavior for an arbitrary po is found to be 

p^^^ ^ p+ + a7]e-'^\ r^oo. (53) 




Thus, the inverse relaxation time is 



and 



2h^/^, h^O, (54) 



h, h CO. (55) 



It is worth making a few comments on the MF approximation. For h = (no input), the 
MF equation predicts in both t ^ decay towards an empty steady state. However, 

if h is finite, the decay becomes exponential with a relaxation time that diverges as /i"^/^ 
for small h. This observation suggests that h = can be viewed as a transition point, 
and one can try to account for spatial fluctuations by some kind of phenomenological 
scaling assumption inspired by the theory of critical phenomena p9|. On the other 



hand, the limit h ^ for fixed R and D considered here automatically realizes the 
continuum approximation a —>■ 0. As expected, ps approaches zero in this limit, but 
the exact solution is larger than MF (pg oc h^^^ vs. pf^^ oc h^^"^). Moreover, the exact 
solution remains larger than the MF approximation over the whole h range ( Fig. 
This is not surprising, since the coalescence step responsible for particle removal requires 
that interacting particles occupy neighboring sites, while in the MF approximation the 
effective range of the interaction is arbitrarily large. On the other hand, as h becomes 
large, the exact steady state approaches the MF curveQ given by Eq. (0) (c.f. Fig. ^. 



*This can also be checked analytically by using the ascending series expansion of the Bessel 
function J^{x) = Efclo iiriJ+k+i) 
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This result is quite natural, since the input mechanism involves no spatial correlations. 
For sufficiently small the system relaxes more slowly into the steady state than the 
MF prediction ( r^^ oc K^^"^ as opposed to t]^^ oc h^^"^ in the MF case). In the large h 
limit, an expansion of the real root of (0) shows that r^j^ oc /i, as expected from the MF 
approximation. The results (|3^ ) for the stationary coverage and (^) for the relaxation 
time in the small input limit are in full agreement with the corresponding continuum 
approximation 0. However, one expects a discrepancy for early times, as observed for 
reversible coalescence. 



III. ANNIHILATION MODELS 

The basic model consists of diffusing point particles that annihilate upon contact. 
The particles hop with probability 2D /a? to nearest neighbor sites; the only difference 
from simple coalescence is that here both particles vanish instantaneously whenever they 
attempt to occupy the same site. This reaction scheme conserves the parity of the 
total number of particles. Parity changes in a given interval may only arise by particle 
migration into or out of the interval. The evolution of Gn{t), the probability that n 
consecutive sites contain an even number of particles, provides detailed information about 
the annihilation process A + A — > S + S. Following similar arguments as those for 



coalescence, one arrives at |j25| 



— — — Gn-i — 2 G„ + Gn+l- (56) 

CLT 

The boundary conditions are Gq = 1 and < < 1. Note the similarity between the 
equations for (for annihilation) and En (for coalescence). Indeed, it has been shown 
that the basic coalescence and annihilation models lie in the same universality class p7 



in the framework of our interval method, one can show that the only difference between 
coalescence, annihilation and the zero-temperature q-state Potts model stems from the 



specific form of the various initial conditions ||25|. In the present case, if the lattice 
has a random distribution of particles with global coverage po, a simple combinatorial 
argument yields 

G„,(0) = ^ + i(l-2po)". (57) 
As predicted by early works of Torney and McConnell [|l^ , and Lushnikov , the 



exact solution of (|56|) (on- and off-lattice) shows an anomalous, asymptotic r~^/^ decay 
of the lattice coverage p = 1 — Gi to the empty steady state. 

The annihilation model may be extended by including input or birth processes. In 
the following, we present exact results concerning some of these possibilities. 
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A. Single particle input 



At each time step a particle is injected at a randomly chosen site at rate Ra. The 
site becomes occupied if it is initially empty and becomes empty otherwise. This is 
equivalent to adding the toggle reactions S ^ A and A ^ S to the original annihilation 
scheme. The relevant set of equations for the probabilities Gn now reads 

= Gn-l -2Gn + Gn+l + U h{l - 2Gn), (58) 

CLT 

with the same boundary conditions as above. Here, the last term represents parity 
changes due to particle input, and r = 2Dt/o? and h = Ra^/2D. 

The model turns out to be closely related to the case of coalescence with input studied 



above. Indeed, Eqs. (|58| ) may be solved with similar techniques to those used in [II B 
The stationary coverage can be expressed in terms of Bessel functions, and for low h it 
is smaller than (|3^) by a factor of 2~^/^. The relaxation time turns out to be smaller 
than that for coalescence by a factor of 2"^/'^ 0. 

These results are in full agreement with previous calculations by Racz, who provided 
an exact solution by mapping the above dynamics into a kinetic Ising model |jl5|. In 
Racz's model, the time evolution of an infinite ID array {cr} = {. . . , cxj, cij+i, . . .} of 
stochastic pseudospin variables ai(t) = ±1 is considered; particles are identified with 
domain walls between regions containing up or down spins only, i.e., the bond variables 
rii = {1 — (jj(Ji+i)/2 are seen as particle occupation numbers. For random homogeneous 
initial conditions, translational invariance holds, and the mean (global) particle coverage 
is given by pit) = (1 — (crjcrj+i))/2. The state probability P({cr}, t) satisfies the Glauber 
master equation 

^^Ml^= £ j:[wi''\{a}t)P{{a}t,t)-wl''\{a})Pi{a},t)], (59) 

i=-oo a=l 

where the state {a}} is obtained from {a} by flipping the ith spin, and {a}j differs from 
{a} by the simultaneous flipping of all spins aj with j < i. The flipping rates are given 
by 

= ^ [1 - ^a,(a.+i + a._i)], wl'\{a}) = Th. (60) 

Setting 5 = 1 and F = 2D /a? corresponds exactly to the dynamics of our model with a 
relative feed rate h. 

It has also been suggested that annihilation with single particle input may be relevant 
to the kinetics of a certain class of processes involving cluster-cluster aggregation in the 
presence of sources and sinks. For example, in aerosol formation aggregation centers 
can be generated by photo-oxidation, and large clusters may disappear as a result of 
sedimentation According to Racz, it should be expected that such models be in the 



same universality class as annihilation with input 
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B. Input of adjacent pairs 



Pairs of particles are injected simultaneously at adjacent sites at rate R per site per 
unit time. Thus, the steps SS —>■ AA, AS — *• SA, SA AS and AA SS take place at 
equal rates (but generically different from the hopping rate for diffusion and annihilation 
events). The kinetics is described by the equations 



dr 



Gn-i ~ 2 Gn + Gn+i + 2/i(l — 2G'„) 



(61) 



with r = 2Dt/a^ and h = Ra'^/{2D). The boundary conditions remain the same as 
before. 

Again, we begin our study of this model with the steady state solution: 



Gr. 



1 1 

2 + 2 



1 + 2/i-a/(1 + 2/i)2-1 



This implies a steady state coverage of 



Ps = yh? + h — h. 
For fast input rates (large h), we have 



1 

8^' 



while for small h. 



Ps 



h oo, 



(62) 



(63) 



(64) 



(65) 



We now derive the exact time dependence of the concentration for a random homogeneous 
initial distribution. Applying the LT to both sides of ( |6TD and using the initial condition 
(P^) yields a set of inhomogeneous difference equations. However, the inhomogeneity does 
not depend on n and can be easily shifted away. After using the boundary conditions, 
we obtain 



Gn 



1+ '1 



2po 



1 



2s 2{s + a-b) \2s 2{s + a - b) 
where a = 2(1 + 2h), b = 1 - 2po + 1/(1 - 2po), and 

s + a- ^{s + aY - A 
A— = . 



A! 



(66) 



(67) 

To obtain the explicit time dependence of the probabilities G„, we invert Eqs. (|66D using 
the convolution theorem for the LT. We thus get 



Gn{r) 



1 n 
2^ 2 Jo 



/n(2r') 



dr' 



r 



:i - 2po)" _n 
2 2 Jo 



-br' 



/n(2r') 



dr' 



■ (68) 
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In the limit po 1/2, the difference equations for the G„ become homogeneous; the 
last term in ( ]68| ) vanishes and one need compute only one remaining integral. 

Determining the long-time asymptotics of the solution ( |68D requires separate analysis 
of two cases: po ^ 1/2 and po > 1/2. As we shall see, both cases yield the same expansion 
for p(r). 

For n = 1, the integrals in the rhs of ( |68|) can be computed by making use of the 
identity |3|] 



-^^/i(2l/) 

y 



oo 



dy 



Z 



2fc+l 



(69) 



where (3; 7) is a generalized hypergeometric function. For short times r, however, 

it is easier to use the series expansion of the Bessel functions to evaluate the integrands. 
In this limit, we obtain, to second order: 

p(r) = Po - 2{pI + h{2po - l))r - {2pl + pl{l + 8h) + Sh^po - ^h^ - + 0(r=^). 

(70) 

The ascending series on the rhs of (|69|) is not suitable for the analysis of the long-time 
asymptotics. Instead, we observe that for po < | {b > 2), the integrals in (|68|) converge 
for r — i> 00. We can therefore split them into a definite and an indefinite part and use 
the adequate asymptotic expansion for the latter. This yields 



p(r ^ 00) = Ps + ^ 



1 1 - 2po 
h pI 



^-3/2 ^-Ahr 



1 



(71) 



For Po > 1 {h < —2), however, the last integral diverges. Nevertheless, we can split the 
integral as follows: 



e'^"'/i(2r') 



dr' 



dr' + r 
Jc 



-br 



'h{2r') 



dr' 



(72) 



/o r' JO r' JCi T 

with a sufficiently large constant Ci, so that the second integral on the rhs can be well 



approximated by expanding the modified Bessel function to its leading term |3^. We 
can then replace the original integral by 



,(2-fe)r' 



1 



+ v/2^erfi(J(2-6)r), 



2^Jc^ r'3/2 

where C2 is a constant and the imaginary error function is defined as 

2 



(73) 



Thus, for n = 1, the last term in 



erfi(a;) = —ieri{ix) 

becomes 



/ dy. 

Jo 



(74) 
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. , D{J{2 - h)T) e^^-'^)- + e^^-"^)- + e^'-'^\ 

2 Jo T V TT ^ 2y/1TT 

(75) 



where C3 = — C2/2 and = e ^ /q^ dy is the so-called Dawson's integral For 
large arguments, 

D{x) = ^ + ^,+o(^). (76) 



e 



(6-a)r 



2x Ax'^ 

If we now use ( |7B| ) to expand the first term on the rhs of ([75|), we end up with 

(77) 

Clearly, the last, purely exponential term decays faster than all other terms, since b < 
—2. Thus, the leading order is the term proportional to r~^/^. We then find the same 
decay law as for po < 1/2, of Eq. ([Til). The relaxation time is universal and given by 

tr = i/m. 

The MF approximation is derived by noting that Gi = 1 — Pi and G2 = 1 — 2 Pi + 2 P2 
and using the factorization ansatz P2 = P^. This leads to: 

^ = _2p2 + 2/i(l-2p) (78) 
dr 

The physical steady state of this equation is given by (|63|) ; this is not surprising, since 
the steady state is an equilibrium one, like in the reversible coalescence case. The early 
time behavior is in agreement with Eq. ( ffOD up to the first order, while the long-time 
approach of the MF solution to ps is easily found to be 

p [T) PsO^^^^^^^^^^y/^^^e , r^oo, [(J) 

Once again, the MF approach fails to capture the time dependent prefactor of the expo- 
nential in the exact solution, although it yields the correct relaxation time for large h. 

In the small h limit, the expressions for ps ~ h^^"^ and tr agree with the continuum 
limit ||2^. These results have also been obtained by Racz by mapping the above model 
into the kinetic Ising model; to do so, one chooses F = 2D jo?, 5 = (1 — /i)/(l -|- h), and 
wf^{\aY) = in Eq. (|59| ) [pTSf . This yields hopping to nearest sites at rate D /a?, nearest- 
neighbor annihilation at rate 2D/(a^(l+/i)), and pair production at rate 2Dh/ {a?{l+h)), 
i.e., the dynamics of our model for small h. At larger /i, however, the correspondence is 
lost because annihilation events in the spin model are slower than hopping onto empty 
sites. Our results are a natural generalization of Racz's solution for the case of arbitrary 
h. In particular, in the large h limit, the steady state coverage is given by Eq. (Bl). 
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We also find that the decay to the steady state is not purely exponential but rather 
given by a power of r times an exponential. This corroborates a recent conjecture of 
Habib et al. |27|, based on the continuum approximation of Eq. (|61|). 

Finally, we remark that the above model might capture essential features of transient 
optical absorption at energies less than the interband gap in certain ID organic semicon- 
ductors |^ , |4T11 . It has been argued that the high energy peak in the absorption spectrum 
of trans-polyacetylene may be due to the photogeneration of intrinsic, self-localized ex- 
citations of the semiconductor chain [4^ in form of alternating soliton-antisoliton pairs 
which move randomly along the chain under the influence of thermal fluctuations and 
annihilate upon contact. 



C. Annihilation with Symmetric Birth 

The input process consists of a double particle birth, i.e., a particle produces off- 
spring at both neighboring sites at rate V (on each side). In other words, we have the 
additional reactions SAS AAA, SAA AAS, AAS -> SAA and AAA SAS. The 



probabilities Gnif) obey the equations P6 



(2D + 2V){Gn+i-2Gn + Gn-i), n>l. (80a) 



dt 

The boundary condition is nontrivial; for n = 1, one has 
dGi 



dt 



2D(G2-2Gi + l) + 2V^(G2-Gi). (80b) 



Notice that in the above equations we have set a = 1. This is done bearing in mind 
that there is no straightforward continuum approximation of the model, due to the form 
of the boundary condition. As before, we have the condition, < G„ < 1. Again, we 
solve the set of Eqs. (|80D for an initially uncorrelated distribution given by the initial 
condition 

We study this model with and without diffusion. Let us first consider the case of 
D = 0. In this limit, clusters of particles can spread only by giving birth. To simplify 
the evolution equations for the Gn we let r = 2Vt and a„(r) = i(l + F„(r)). We then 
have the initial BVP: 

dF 

= Fn+i - 2Fn + Fn-1, (81a) 

dr 

-7^ = ^2 - Fi, (81b) 
or 

F„(0) = (1 - 2po)". (81c) 

The solution of these equations is readily obtained by LT methods similar to those used 
above. The LT of F\{f) takes the form 
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Flis) 



l-2po 



Po 



Is + A 



V 



V 



12) 



where v = 4pq/(1 — 2po)- For three simple cases, po = 0, po = 1/2, and po = 1, the 
transform can be inverted to yield Gi = 1, Gi = 1/2, and Gi = 1/2(1 — e^^'^/o(2r)), 
respectively. These results merely indicate that an empty system remains empty; po = 
1/2 is the steady state coverage; and initially full lattices decay into the steady state. 
We can also invert explicitly Fi{s) when < po < 1/2. We get 



Fiir) 



Po - Po{v + 



4) re-("+2)r'j^(2r')rfr' 
Jo 



Poe-^Vo(2r) 



Using the appropriate long time expansion, we get 

l-2po 



Fir) 



PoV4 



r — > oo, 



TCT 



or, in terms of the coverage 
Pir) 



1 1 J, , , 1 1 - 2po 

Fiir) ^ r=, T 

2 2 ^ ^ 2 4pov^ 



oo. 



^3) 



^5) 



Thus the system decays algebraically into the steady state p^ = 1/2. In fact, Sudbury had 
shown that any homogeneous, random distribution leads to a half-empty lattice |l30| , p6| . 
It is therefore not surprising that we obtain the same behavior also for 1/2 < po < 1. 
In this regime, t; < 0, and the LT inversion becomes difficult. Since we are mainly 
interested in the long time asymptotics of Fi, we circumvent this difficulty by expanding 
Fi{s) about s = and making use of a Tauberian theorem. For s — > 0, 



Fi{s) 



2po 



implying that 



FAr) 



l-2po 
PoV4 



r 



oo. 



ITT 



This matches the result for < po < 1/2. 

We now turn to the case of annihilation with double birth and diffusion. 
Again, we take Eqs. ( ^OD and (pTj) as a starting point; setting r = 2{D + 
Gn = ^{Fn + 1), and taking the LT with respect to r, we finally obtain 



^6) 



(87) 



D > 0. 
V)t and 




2po(l - k) + 



+ As-s-2k\ 1 - 2po 



s — V 



(1 — k)s — 



S — V 



with K = D/[D + V). Notice that we recover the appropriate solutions when k = 1 and 
K = 0. These limits refer to the V = and D = cases, respectively. Fi{s) can be 
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inverted for particular cases of the parameters po and k. However, it is simpler to make 
use of the s = expansion to obtain the long time asymptotics: 

s^O. (89) 



This leads to 



This result is not valid for k = 0, since the limit k — > is singular. It is easy to see that 
the system cannot become empty without diffusion pO| , ^ . In that case, ( pSf ) holds. 



p^-^, r^oo. (90) 



Let us compare the above results with the MF approximation, which can be derived 
from the equation for Gi in the usual way. In the presence of diffusion, the simple MF 
approach yields 

g = (1 - «:)p - 2p^ (91) 

Consequently, the steady states are ps = 0, (1 — k)/2. Note that the existence of a non- 
empty steady state is already in contradiction to the exact solution (^). The nontrivial, 
time dependent solution of (|91|) reads 



P 2po + (l-«:-2po)e(-i)^' ^^^^ 

predicting exponential relaxation, in contrast to the inverse power law behavior of the 
exact solution. Thus, for k > 0, the MF approximation fails both at the static and the 
dynamic levels, whereas for k = (no diffusion), it describes the steady state correctly, 
but not the long time dynamics. At the MF level, the effect of diffusion consists in 
driving the system from a steady state with half the lattice filled into a steady state with 
a smaller (but finite!) coverage. In contrast, the effect of diffusion on the exact solution 
is much more drastic, since the slightest amount of mobility already lands the system in 
an empty absorbing state. Note also the somewhat surprising fact that the presence of 
diffusion drives the steady state away from the MF prediction, except when the value of 
K becomes close to one. In this limit, p^^^(oo) again becomes arbitrarily close to zero. 



IV. SUMMARY AND OUTLOOK 

We have derived a series of exact results for the ID kinetics of the lattice coverage and 
the particle concentration for coalescence and annihilation with various particle sources. 

In the case of reversible coalescence the MF approach reproduces the equilibrium 
steady state correctly, but it only yields the relaxation time of the exact solution if h and 
Po are large enough. Besides, this approximation is not sensitive enough to reproduce the 
observed fluctuation- induced phase transition. In contrast, the off-lattice approximation 
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valid for low input rates h (and thus for low transient concentrations) provides a good 
qualitative description of the long time asymptotics. For large /i, no qualitative changes 
in the cq and r dependence are observed, but the values of the prefactor and the relaxation 
time in the asymptotic form of the solution depend on the lattice constant a. For early 
times, we find that the exact solution is in agreement with the classical MF approach, 
while the off-lattice approximation dictates a slower kinetics in this regime (one expects 
a similar conclusion for the other lattice reactions studied in this paper). 

In the case of coalescence and annihilation with single-particle input, the exact steady 
state is expressible in terms of Bessel functions. For small /i-values, it is well reproduced 
by the continuum approach, whereas for large values, it becomes MF-like. In this case, 
the relaxation is purely exponential, and tr is again MF-like only for large h. 

For annihilation with pair input, the steady state is an equilibrium state and can be 
derived directly from the MF approach. However, the relaxation is not purely exponen- 
tial, as predicted by the MF approach, since the concentration behaves like r~^/^e~^'''^^ 
with a universal tr. In this case, the somewhat faster relaxation with respect to the 
single particle case may be due to the higher particle input. 

Summarizing, in the cases for which the continuum approximation was considered, 
the latter is in good qualitative agreement with the exact long time kinetics and with 
the steady state in the low h limit, however, deviations from the exact kinetics are 
expected for large values of /i, i.e., at higher transient concentrations. In this limit, the 
MF approximation yields relatively good results for even in the reversible coalescence 
case, provided that the initial concentration is sufficiently high. Fluctuations can also 
be safely neglected in the short time regime. 

As for annihilation with symmetric birth, we were able to confirm Sudbury's result 
for an infinite lattice, and prove that the long time relaxation to the steady state is 
proportional to r~^/^, both for immobile reactants and in the presence of diffusion. In 
the diffusionless case, the MF approach reproduces the exact steady state correctly, but 
not the long time kinetics. Any amount of diffusion lands the system in an empty state, 
in contrast to the MF prediction. 

Possible generalizations of our work include studying finite size effects and general- 
izing the method to better understand the spatial organization in these systems; in this 
respect, it is of interest to compute quantities such as interparticle distribution func- 
tions or multiple-point correlations from joint probabilities for pairs of intervals, both 
on-lattice and in the continuum limit WM- The present methods also allow to study the 



effect of temporal correlations |20] or inhomogeneities in the initial conditions 
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FIGURES 




FIG. 1. Exact lattice coverage (solid lines) vs. mean field approximation (dashed lines) of 
the reversible coalescence for two different initial conditions po = 0.75 and po = 0.2 {h = 1). 
Since 0.2 < 0.29 ~ pc, the exact solution relaxes much slowlier into the equilibrium steady 
state in the latter case 




FIG. 2. Comparison of the approximations of Ei obtained by truncating the series expan- 
sion (|40|) to 1st (dashed), 2nd (dashed-dotted) and 3rd order (solid) with the exact solution 
(diamonds). The relative feed rate is h = 1. 
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FIG. 3. Comparison of the steady state coverage as a function of the relative feed rate h 
computed from the exact solution and the MF approach. 
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